Dynamical phase transition of a periodically driven DNA 
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Replication and transcription are two important processes in living systems. To execute such 
processes, various proteins work far away from equilibrium in a staggered way. Motivated by this, 
aspects of hysteresis during unzipping of DNA under a periodic drive are studied. A steady state 
phase diagram of a driven DNA is proposed which is experimentally verifiable. As a two state 
system, we also compare the results of DNA with that of an Ising magnet under an asymmetrical 
variation of magnetic field. 



I. INTRODUCTION 

The semi-conservative replication of DNA requires 
a complete separation of the two strands, which oc- 
curs concomitantly with the production of new strands. 
Even transcription needs a partial separation of the two 
strands of DNA to read the genetic code of the base se- 
quence [l| . The traditional view of a temperature or pH 
induced melting is getting superseded by a more mechan- 
ical procedure of force or torque induced unzipping 0- 
[sj . This paradigm shift is because of the lack of extreme 
temperatures or the environment required in general for 
the melting (the melting temperature is around 70-lOOC) 
and the recent establishment of the unzipping as a gen- 
uine phase transition [Ij. A double-stranded (ds) DNA re- 
mains in the bound doublet state below its melting tem- 
perature even when pulled by a force at one end until and 
unless the force exceeds a critical force g — gc{T), above 
which it opens into two single stranded (ss) DNA. Many 
aspects of this transition have been studied since then 
both theoretically [^-[l^ and experimentally [1, 16l [it} . 
mostly in equilibrium, though puzzles remain |18. Il9|. 

In vivo, a DNA is opened by helicases, which are mo- 
tors that move along the DNA[H. Both the motion and 
the opening processes require constant supply of energy. 
A few examples are DNA-B, a ring like hexameric he- 
licase that pushes through the DNA like a wedge [20|. 
PcrA that goes through cycles of pulling the ds part of 
the DNA and then moving on the ss part [2l| , viral RNA 
helicase NPH-II that hop s cy clically from the ds to the 
ss part of DNA and back|22|. Such cycles of action and 
rest, with the periodic ATP consumption, indicates an 
exertion of a periodic force on the DNA. 

It suffices to describe the equilibrium unzipping transi- 
tion by the two thermodynamic conjugate variables, force 
g and separation x of the pulled base pair (Fig.l). In 
thermal equilibrium, a quasi-static change of the force 
from zero to a maximum gm and then back to zero, keep- 
ing other intensive quantities fixed, would result in re- 



* |garim abhu@gmaiI. com| 
t |poulomi@iopb. res.in| 

somen@iopb.r es.in| 
S lyashankit @yahoo . com | 



tracing the thermodynamic path, ending at the initial 
state; history plays no role because of ergodicity. In con- 
trast, under a periodic force, the mismatch between the 
relaxation time and the external timescale for change of 
force would restrict the DNA to explore a smaller re- 
gion of the phase space, creating a difference in the re- 
sponse to an increasing or a decreasing force. Near a 
phase transition, where the typical time scales of dynam- 
ics become large, the difference between the forward and 
the backward branches becomes prominent. This is hys- 
teresis of DNA unzipping (23|. More recently Kaprip^l 
showed how the work theorem can be used via a multi- 
histogram method to extract the equilibrium isotherm 
from the hysteresis curves. The nature of hysteresis and 
its dependence on the applied field and frequency are 
well studied in the context of magnetic and structural 
systems [2^, [2^, though it is the time- averaged loop that 
has received the attention. With the advent of single- 
molecule experiments on short DNA chains (oligomers), 
it might be possible to probe the time-resolved loops, not 
just averages. Motivated by the biological relevance and 
the experimental feasibility, we explore the behavior of 
DNA under a periodic force, to be called a periodically 
driven DNA. Our results show that without changing the 
physiological conditions {e.g. the temperature or pH of 
the solvent), a DNA chain may be brought from the un- 
zipped state to the zipped state and vice versa by varying 
the frequency (i^) alone. By using a probabilistic descrip- 
tion of the time variation of the DNA response, we pro- 
pose a force-frequency (g-iy) phase diagram for the driven 
DNA. 

A similar dynamic symmetry breaking transition is 
known to occur in ma gne ts where the nature of the 
hysteresis loop changes [23. Because of the two sta- 
ble phases, the unzipping transition has often been de- 
scribed by a two state Landau type free energy func- 
tional 0, [231 . As a two phase system, we make use of 
the magnetic Ising model to corroborate the behavior of 
DNA, both undergoing a first order transition and show- 
ing hysteresis. We establish that the observed features 
and the phase diagram are robust and generic, and not 
tied to any particular DNA model. 

It was recently shown that the DNA hairpin (Fig. 1) of 
32 beads (stem 10 base pairs and loop of 12 bases) cap- 
tures the essential properties of a long DNA chain with 
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FIG. 1. DNA hairpin in a zipped (Z) and an unzipped 
(U) state. The stem (soUd lines) consists of complementary 
nucleotides, whereas the loop (dashed line) is made up of non- 
complementary nucleotides. Only native interaction (dotted 
lines) for base pairing is allowed. For this paper, the stem is 
of length p — 10, and the loop has 12 monomers. The total 
length is P = 32. 



middle of the chain is experimentally difficult. More- 
over, the biological situation requires the force following 
a moving replication fork, vifhich is also not amenable to 
single molecule force experiments. In view of these, we 
use a simpler geometry where the force is applied at one 
end of the liairpin keeping the other end fixed, which 
could be used in single molecule experiments. The typi- 
cal time scale of fluctuations involved in unzipping varies 
from ns to /zs, and therefore we avoid the details of an all 
atom simulation by adopting a minimal ofF-lattice coarse 
grained model of a DNA hairpin, where each bead repre- 
sents nucleotide as a basic unit, which comprises a base, 
a sugar and a phosphate group. We have performed our 
simulation in reduced unit as discussed below. 



implicit bubbles [23| . The force-temperature diagram of 
the DNA hairpin is found to be qualitatively similar to 
the DNA unzipping experiment |16| in the entire range of 
/ and T, whereas the phase diagram for a dsDNA of 16 
base pairs differ significantly with the experiment. Since, 
under the periodic force, the chain may unzip and rczip, 
the bubble dynamics and the wandering interface (Y-fork 
junction) meeting transient bubbles (extra source of en- 
tropy), would play important roles in hysteresis. Need- 
less to say that such bubbles are ubiquitous in a long 
chain. A DNA hairpin is the simplest model, where the 
premadc loop allows us to study the dynamics of such 
implicit bubbles. 

The outline of the paper is as follows. In Sec. |lll we 
introduce the model used for the hairpin and the method 
of study. The results arc analyzed in terms of average 
hysteresis loops and time resolved loops. The dynamic 
phase diagram can be found in this section. In Sec lIIIl the 
hysteresis loops of an Ising ferromagnet under an asym- 
metric field modulation is given. The similarities with 
the DNA hairpin problem are discussed here. The paper 
ends with a short summary and conclusion in Sec. IIVI 



II. DNA HAIRPIN UNDER PERIODIC FORCE 

A. Model and method 

DNA hairpins (see FigH]) occur naturally in vivo and 
are used often in experiments on DNA with bubbles [l^ . 
The non-paired bases of the loop is a source of entropy 
which in turn controls the dynamics of hysteresis. In 
this paper, we perform Langcvin dynamics (LD) simula- 
tions of a DNA hairpinflj, to study the separation x of 
the terminal base pairs pulled by a periodic force g(t) of 
time period t(= l/iy) (Fig. [5]). There have been efforts 
[H, [l^ to understand transcription like processes by ap- 
plying a constant force in the middle of a DNA chain. 
A very rich phase diagram has been obtained though no 
such experiment has been reported so far to confirm the 
phase diagram. It appears that applying a force in the 



1. Model 

The configurational energy of the system under con- 
sideration is [2 3,] 

p p 2 P 

= 5: - rfof + E E 4(;^ -#)•(!) 

4 = 1 1=1 J>i+1 »J 

Here, P = 32 is the number of beads. The harmonic 
term with a spring constant k (=100) couples the adja- 
cent beads along the chain. The distance between the 
adjacent beads dij's defined as jr^ — rj], where and rj 
denote the positions of beads i and j, respectively. We 
assign C = 1 of the Lennard- Jones potential. The base 
pairing interaction Aij = 1 is restricted to native con- 
tacts only, which is similar to the Go model. All pairs 
of nucleotides that do not form native contacts in the 
stem and in the loop interact only through short range 
repulsion (excluded volume). The parameter do(= 1.12) 
corresponds to the equilibrium distance in the harmonic 
potential, which is close to the equilibrium position of the 
average L-J potential. We obtained the dynamics of the 
system by using the following Langevin equation [2^ [s^l 

cPr dr „ ^ 

where m is the mass of a bead which is set equal to 
one here, and the friction coefficient used in simulation 
is C = 0.4. Fc is defined as the gradient of the energy, 
— ViJ, and r is a random force, a white noise with zero 
mean and correlation 

< T,{t)T,it') >= 2CkBT5,j{t - t'), (3) 

where fc^ is the Boltzmann constant which is set to one 
in present simulation. The choice of Eq. |3] keeps the 
temperature of the system constant through out the sim- 
ulation for a given force. In order to study the behaviour 
of a DNA hairpin under a periodic force, we add an addi- 
tional energy —g{t).x{t) to the total energy of the system 
EqlT] We have performed the simulation at T = 0.10, for 
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50 100 150 200 

No. of LD step 

FIG. 2. A typical periodic force used in the simulation. A 
cycle refers to one period of variation of the applied force g 
from to (called amplitude) and back to at a fixed time 
period t. This example has amplitude Qm = 0.1 and two 
dilTerent frequencies u — 1/t — 0.38 and 0.19 for n=5 and 10 
respectively. 



2. Integration 

The equation of motion is integrated using a sixth 
order predictor corrector algorithm with a time step 
A<=0.025 llOl- Ref. [l^l showed that the present coarse- 
grained model exhibits the observed features of the ex- 
perimental g-T diagram [l^ of a long DNA with implicit 
bubbles as well as hysteresis at low temperatures in the 
g-x plane as seen experimentally [l7j . During the simu- 
lation, g is changed in steps of Ag(= 0.01) from to a 
maximum and then to (Fig. This one period 
is to be referred to as a cycle and gm as the amplitude. 
The time period 



9n 



T = 2nAt(l + ^), 
Ag" 



which gc ^ 0.20 [3l|. Wc may convert the dimcnsionless 
units to real units by using the relations: (29j. 



t* = 



1/2 



t, r 



(4) 



is controlled by n, the number of LD steps executed after 
every increment (or decrement) of force. We choose n 
such that nAt is much below the equilibration time. 

By changing g„i or i^, we find it is possible to induce 
a dynamical transition between a state of time averaged 
zipped (Z) or unzipped (U) to a dynamical state (D) 
oscillating between Z and U. 



where t*, r* and e are time, distance and characteristic 
hydrogen bond energy in real units, tr is the distance 
at which inter-particle potential goes to zero. For exam- 
ple, if we set effective base pairing energy e ~ 0.1 eV, 
the average mass of each bead 5 x 10~^^ g and ct = 5.17 
A, we get a time unit r ss 3 ps. It is known that ^ 
affects the kinetics only and thermodynamics remain un- 
changed [13, [s^l ■ The friction coefficient used in our sim- 
ulations C = 0.4mT~^ = 6 X lO^^^gs^^. The order of 
force thus would be e/A~ 160 pN. For temperature con- 
version, a little more care is needed because the coarse 
grained model does not take into account many other 
energies 3^ . We compare the phase diagrams for DNA in 



the thermodynamic limit using a two state model based 
on the modified Freely Jointed Chain (mFJC) model of 
polymer with the experiment [l6j . The melting temper- 
ature at 363K and the flat portion (the stretching of co- 
valent bonds) at 293K in experiments [l^ correspond to 
0.23 and 0.05 in reduced units, which are also consistent 
with the simulation (23| . If we use this information to 
scale the temperature in a linear fashion, we have a rela- 
tion, T* = 363 -h^||5§^(r- 0.23), where T* is the real 
temperature. This mapping is valid only in the range of 
overlap of mFJC and experiment. This relation gives the 
reduced temperature used here T = 0.1 as equivalent to 
a real temperature around 320K, i.e., around 45C. It may 
be noted that h yst eresis has been observed in this range 
of temperature [l7| . 



B. Numerical analysis 



Hysteresis loops 



In Fig 131 we plot the average value of x{g) over C (= 
1000) cycles vs 5, for different values of gm and v. These 
loops for different initial conformations remain almost 
the same (except (23)) so that the large menagerie of 
shapes of hysteresis loops observed are typical, not ex- 
ceptional or accidental. At a high frequency, the DNA 
remains in Z or in U depending on whether gm < '2gc 
or not, (Figl3J13)&(43)), irrespective of the initial con- 
formation. For a relatively smaller j/, the sequence of 
hysteresis loops (FiglSl (11,21 & 31)) is of a different na- 
ture. These clearly reflect that the DNA starts from the 
Z (x = at the start of the cycle), goes to the U (x = 30) 
and then back to the Z [sl]. Some of these plots show 
the phase lag between the force and the extension, e.g., 
even when the applied force decreases from gm to (Fig. 
131 (32 & 33)), the extension x{g) increases. If the sys- 
tem gets enough time to approach equilibrium, the lag 
disappears (e.g., Figl3][21 & 31)). A different behavior is 
seen in the case of intermediate forces (Figl31 (21, 22 & 
23) ). Despite the identification of the states at the two 
extreme forces as Z and U, there is a significant change in 
the X- values ed, g = and g = gm, in Figl31j23), compared 
to the other loops shown in Figl^l Most striking here is 
the wide sample to sample fluctuations in the loop. 
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FIG. 3. (Color online) DNA hysteresis for different gra and v (increases from left to right) as indicated. Each plot contains 
the loops for 10 different initial conformations. These are at T = 0.10, for which gc ^ 0.20. 



2. Time-resolved loops 

The failure of the average response to provide a de- 
scription of the steady state dynamic behavior prompts 
us to analyze the distribution of paths over the different 



cycles in terms of a dynamical order parameter Q defined 
as 



Q = - 

T 



x(t) dt. 



(5) 
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This Q is analogous to the dynamical order parameter 
studied in the context of magnetic systems [26j . In FigHJ 
we plot the values of Q over different cycles for three 
different v. The time sequence of Q seems not to indicate 
any regular pattern, and therefore, we assume that the 
allowed states occur randomly. The time sequence can 
then be interpreted in terms of a probability of getting 
a particular value of Q. We find that the steady state is 
described by a stationary probability distribution {P{Q)) 
which are also shown in Fig The plots show three 
peaks which could be associated with phases U, D and Z 
by comparing with Fig. [31 Based on the width oi P{Q) 
we make operational criteria, < Qz < 10, Qu > 22, 
and 10 < Qd < 22 for the identification of the phases. 

From the peak locations of P((5), we map out the phase 
diagram of the driven DNA in the g-iy plane (Fig [5]). 
A line in FiglS] represents a boundary beyond which a 
particular peak appears or disappears and resembles a 
first-order line. These are now discussed sequentially. 

(A) For a very small force, P{Q) expectedly shows 
a peak in the Qz window for any v (Fig. |4] 1st row). 
Beyond a certain gm there appears a second peak in the 
Qd window. See Fig. H 12 to 22 or 13 to 23. For each 
chosen v, we determined the value of gm above which the 
peak of Qd appears, giving us the lower boundary for D. 

(B) For a given v, with increasing g^, the peak height 
of Qz decreases as seen in Fig. |4] along any vertical se- 
quence, top to bottom. One gets one or two other peaks 
but Z disappears beyond a certain gm, giving an upper 
boundary for Z. 

(C) Before Z disappears, there could be peaks for D or 
U, or both ( Fig. 111(31, 22)), implying coexistence, Z+\J 
or Z-l-D-l-U, for a range of gm- From the appearance of 
U, we determine the lower boundary for U. 

(D) Then, as Z disappears, the system may be either 
in the D phase (Fig. |4l 23R), or, in the mixed state of 
D+V (Fig. m 32R). To be noted is the crossing of the 
boundaries of U and Z, that produces the region of D 
phase. Once the mixed phase D-l-U appears, a further 
increase of gm would vanish the D peak and only U peak 
survives. 

An important point is the occurrence of different Q val- 
ues or regions of phase coexistence. In such a situation, 
the average response is not a good measure of the state 
and there will be intrinsic difficulty in reducing fluctua- 
tions in the hysteresis curve. 



III. DYNAMICAL PHASE TRANSITION IN 
ISING FERROMAGNET 

To use the magnetic analogy, consider a ferromagnet 
like a two dimensional Ising model below its critical point, 
subjected to a periodic magnetic field h = ho cos(2'Kh't) 
in time (t). As v increases, there is a dynamic transition 
where the average magnetization in a cycle goes from 
zero to a nonzero value p6j. For small v there is a hys- 
teresis loop connecting the two symmetrically opposite 



magnetized states, while for large v, the lethargic sys- 
tem remains magnetized in one direction. The analogy is 
between the conjugate pair [h, m), where m is the magne- 
tization, for a magnet and {g, x) for DNA. The magnetic 
hysteresis loop in the h-m plane is analogous to that in 
the g-x plane. To simulate a DNA-like behavior in the 
Ising system, an asymmetric modulation of the magnetic 
field has been taken to make the average field over a cycle 
different from the critical value. A Monte Carlo dynam- 
ics is used to study the hysteresis in an 8 x 8 square lattice 
Ising model with periodic boundary conditions under a 
periodic field between hi and /i,„. 

A. Model and Method 

We consider a two-dimensional Ising spin system 
({s,; = ±1}) with nearest neighbor interaction and un- 
der a magnetic field /i, 

H =-JY,S^S,-hY,S^, (J>0), (6) 

with i labeling the sites of an 8 x 8 square lattice with 
periodic boundary conditions. The infinite lattice crit- 
ical temperature is Tc « 2.269J/fcB in zero field. The 
magnetization is defined as the thermal average 

m = N-'Y.^s,)^ (7) 

i 

of N{= 64) spins. For the above Hamiltonian, we choose 
J = 1 and ksT = 2 with fcs = 1, so that T < T^. 

The Monte Carlo procedure used is as follows. We 
choose a spin, calculate the change in energy AE of 
the system if the spin gets flipped. Whether this 
spin would be flipped or not is determined by us- 
ing the standard Metropolis algorithm, which gives 
the probability of acceptance of an attempted flip by 
^accept = niin(l, e^'^'^^). In this way, we sequentially 
consider all the N spins, one at a time, to attempt a 
flip. The time taken to attempt N spin flips constitutes 
one MC sweep. 

A complete cycle consists of 2Af steps, starting from 
hi, reaching hm and then back to hi. Initially system 
is equilibrated at hi = —0.6 and ksT = 2. Then the 
periodic magnetic field is switched on. At each step, (i) 
the magnetic field is increased by Sh = {hm. — hi)/M and 
(ii) the system is relaxed towards equilibrium by n = 5 
MC sweeps, which is much less than the equilibration 
time. The magnetization m is calculated at each of such 
27V steps. The average of magnetizations calculated over 
a cycle then gives the quantity Q. The above process is 
repeated C = lO'' times, i.e., for 10^ cycles. The Q values 
obtained in this way gives the probability distribution 
P{Q) for a given hm- We simulate the system for various 
frequencies (controlling M) and hm, keeping hi fixed. 

An asymmetric field in the Ising model enables us to 
distinguish two differently ordered phases, the counter- 
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FIG. 4. The time sequence of Q and P{Q) vs Q plots. Three pairs of vertical panels for three i^'s. Each panel has (i) the time 
sequence on the left (L), and (ii) P(Q)-plot on the right (R) for different (?m's. 



parts of the zipped and the unzipped states, and in addi- 
tion a hysteretic state, to be called the dynamic state D. 
For easy comparison, the negatively magnetized state is 
named Z while the positively magnetized state is U. The 
states are distinguished by a dynamic parameter, 

Q = - f m(t) dt, (8) 

T Jo 

and the operational definition adopted is (i) Qz = 
{-64 < Q < -40}, (ii) Qu = {55 < Q < 64}, and 
(iii) the rest is Qd- Cases (i) and (ii) occur when the 
paths in the m-h diagram remain on one side throughout 
the cycle, and (iii) the dynamical state, D, occurs either 
if there are paths connecting positive and negative values 
of magnetization or if the paths remain more or less near 
zero of magnetization. The division of the three intervals 
or regions are independent of hm ■ 
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FIG. 5. (Color online) Dynamical phase diagram of a peri- 
odically driven DNA hairpin in the Qm-i' plane. The lines 
are boundaries for various phases U, D and Z. The points are 
from the simulation and the lines are guide to eye. 



B. Numerical results 

The results for P{Q) vs Q shown in Fig. |6]are very sim- 
ilar to Fig. Uland the interpretations toe the DNA line. 
We point out a few salient features here. For /im < 0, 
there is only the Z state (Fig. |B^bc). With increasing 
hrrn P{Q) dcvclops a multi-pcakcd structure which indi- 
cates coexistence. Various sequences of states seen are 
noted in the figure caption. These include D-l-U, Z+D 
and Z-l-D-l-U. We find a clear dependence on v of the oc- 
currence of the various states, similar to the DNA case. 
There are situations where U appears before Z vanishes. 
In this regime, we see the coexistence of three phases, 
Z-HD-hU (e.g., = 0.6 in FigE]). The peak height for D 
first grows from zero and then depending on frequency, 
it may decrease, with a gradual appearance of a U peak 
(Fig. |6^,b), or it may keep on increasing up to a very 
high field seemingly merging with U. (low frequency case 
of Fig. int). The Ising case has a special line at hm = —hi 
where the + ^ — symmetry is restored. On that line, 
there is the dynamic transition (D-f->- Z or U) at a partic- 
ular frequency [2^. The occurrence of the Z-l-D-l-U coex- 
istence region found on this line is a finite size effect and 
is a novelty of the mesoscopic system. The h-iy phase dia- 
gram is qualitatively similar to that of the DNA hairpin, 
shown in Fig [51 



IV. SUMMARY 

This paper reports the possibility of a dynamical tran- 
sition of a DNA hairpin and a magnetic spin system both 
under a periodic drive. When a DNA is subjected to 
a periodic force, or, a magnetic system is subjected to 
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FIG. 6. (Color online) P{Q) vs Q for the Ising model for three different v's, and the same set of hm with hi = —0.6. For a 
given v, the distribution P{Q) shifts as hm increases. It goes from Z to U through D. The sequence of phases are (a) Z — >■ Z+D 
Z+D+U ^ D+U ^ U. (b)Z ^ Z+D ^ D ^ D+U. (c) Phases: Z ^ Z+D ^ D. 



a periodic magnetic field, then the system driven away 
from equilibrium shows hysteresis. Usually to get the 
hysteresis curve, one averages over many cycles of the 
field or force. But we show here that this averaging sup- 
presses or loses the actual picture of the states. This 
can be observed if one looks at the individual cycles. By 
quantifying the state of a cycle by Q, the average sep- 
aration for DNA or the average magnetization for the 
spins, we find that the distribution of Q is multi-peaked. 
Whenever a distribution is multi-peaked, the average is 
not a good representation of the state. We identify the 
peaks of the distribution with the zipped (Z), unzipped 
(U) and the dynamical (D) state for the DNA. There is 



a dynamical transition from one phase or state to an- 
other phase, and also there are coexistence of phases in 
some range of strength of drive and frequency. The phase 
diagrams constructed from the data for both the DNA 
hairpin and the magnetic system are of similar nature. 
The role of hysteresis in biological processes remains an 
unexplored territory. On a more quantitative front, our 
work calls for time resolved experiments on periodically 
driven DNA hairpins to explore such hitherto unknown 
dynamical phase diagrams. 
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